library(SDMTools)
wd='/home/jc148322/rob/inputs/';setwd(wd)

#read in the asciis & create base.asc
bdall=read.asc('Bdall.asc')
bdwet=read.asc('Bdwet.asc')
bc20=read.asc('bc20.asc')
canopy=read.asc('canopy.asc')


#read in the data
pos=read.csv('canopy_bc20.csv',as.is=TRUE)
threshold=read.csv('maxentResults.csv')
threshold=threshold$Minimum.training.presence.logistic.threshold
pos.frogs=read.csv('frog_pos.csv', as.is=TRUE)

##################################################################################################################
#CREATE THE LAYERS, EDIT THE DATA
out.dir='/home/jc148322/rob/outputs';setwd(out.dir)
#create the sig.diff ascii
bdall[which(bdall<threshold[1])]=0; #set values below threshold to 0
bdwet[which(bdwet<threshold[2])]=0; #set values below threshold to 0
######################################
#subset frog occurrence data and cbind bc01 and bc12
pos.frogs$bc20= extract.data(cbind(pos.frogs$Easting,pos.frogs$Northing),bc20)
pos.frogs$canopy = extract.data(cbind(pos.frogs$Easting,pos.frogs$Northing),canopy)

frogdata=pos.frogs[,c('Species','bc20','canopy')]
lorica=frogdata[which(frogdata[,1]=='LITLORI'),];
nannotis=frogdata[which(frogdata[,1]=='LITNANN'),]; nannotis=nannotis[,2:3]

record.time=c('New','New','New','New','New','New','Old','Old','Old','Old','Old','Old','Old','Old','Old')
lorica=cbind(lorica,record.time)

######################################
#edit and simplify the data for biplot
tdata = cbind(pos$bc20,pos$canopy); #

tdata = unique(tdata) #find only unique combinations of rounded bc20 and canopy

pos$bdall[which(pos$bdall<threshold[1])]=0; #set values below threshold to 0
pos$bdwet[which(pos$bdwet<threshold[2])]=0;

wet.points=cbind(pos$bc20[which(pos$bdwet>0)],pos$canopy[which(pos$bdwet>0)])
all.points=cbind(pos$bc20[which(pos$bdall>0)],pos$canopy[which(pos$bdall>0)])
		
##################################################################################################################
#CREATE PLOTS AND IMAGES
plotcols=c('grey','tan','forestgreen')
#make the biplot
png('biplot_canopy_bc20.png', width=8,height=8,units='cm', res=300,pointsize=5)

        plot(tdata, xlab="Mean Annual Solar Radiation",ylab='Canopy cover %',xlim=range(tdata[,1],na.rm=T),ylim=range(tdata[,2],na.rm=T), type='n',cex.lab=1, cex.axis=1, font.lab='2')
        points(tdata, col='grey', pch=16,cex=1.5)

        points(all.points,col='tan', pch=16,cex=1.5)

        points(wet.points,col='forestgreen',pch=16,cex=1.5)

        points(nannotis,col='grey20',pch=4,cex=0.7,lwd=0.5)
				
        points(lorica$bc20[which(lorica$record.time=='New')],lorica$canopy[which(lorica$record.time=='New')],col='black',bg='red',lwd=0.5, pch=21)
		points(lorica$bc20[which(lorica$record.time=='Old')],lorica$canopy[which(lorica$record.time=='Old')],col='black',bg='orange',lwd=0.5,pch=21)
		
		legend(23,7000, c('Climate space', 'New Bd SDM', 'Original Bd SDM'), fill=plotcols, cex=1, bty='n')
		legend(23.1,6000, c('L. nannotis', 'L. lorica (old records)', 'L. lorica (new records)'), col=c('gray20','black','black'),pt.bg=c('red','orange'),pch=c(4,21,21), pt.lwd=0.5, cex=1, bty='n')
		
dev.off()

#nannotis points within wet and all
		poly.wet=chull(pos$bc20[which(pos$bdwet>0)],pos$canopy[which(pos$bdwet>0)])
		poly.wet=c(poly.wet,poly.wet[1])
		poly.all=chull(pos$bc20[which(pos$bdall>0)],pos$canopy[which(pos$bdall>0)])
		poly.all=c(poly.all,poly.all[1])
		
		nann.wet=pnt.in.poly(nannotis,wet.points[poly.wet,])
		nann.all=pnt.in.poly(nannotis,all.points[poly.all,])
		nann.dry=sum(nann.all$pip)-sum(nann.wet$pip)
nann.count=as.data.frame(rbind(c('nannotis.all','nannotis.wet','nannotis.dry'),c(sum(nann.all$pip),sum(nann.wet$pip),nann.dry)))
write.csv(nann.count,'nannotis.count.canopy.csv', row.names=FALSE)
